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A novel approach to the description of superconductors in thermal equilibrium is developed within 
a formally exact density-functional framework. The theory is formulated in terms of three "den- 
sities": the ordinary electron density, the superconducting order parameter, and the diagonal of 
the nuclear TV-body density matrix. The electron density and the order parameter are determined 
by Kohn-Sham equations that resemble the Bogoliubov-de Gennes equations. The nuclear density 
matrix follows from a Schrodinger equation with an effective TV-body interaction. These equations 
are coupled to each other via exchange-correlation potentials which are universal functionals of the 
three densities. Approximations of these exchange-correlation functionals are derived using the di- 
agrammatic techniques of many-body perturbation theory. The bare Coulomb repulsion between 
the electrons and the electron-phonon interaction enter this perturbative treatment on the same 
footing. In this way, a truly ab-initio description is achieved which does not contain any empirical 
parameters. 

PACS numbers: 74.25.Jb, 74.25. Kc, 74.20.-z, 74.70. Ad, 71.15.Mb 



I. INTRODUCTION 

One of the great challenges of modern condensed- 
matter theory is the prediction of material specific prop- 
erties of superconductors, such as the critical tempera- 
ture T c or the gap at zero temperature Ao. The model 
of Bardeen, Cooper and Scrieffer (BCS)i successfully 
describes the universal features of superconductors, i.e. 
those features which all (conventional, weak-coupling) su- 
perconductors have in common, like the universal value 
of the ratio 2Ao/fcBT c . The great achievement of BCS 
theory was the microscopic identification of the super- 
conducting order parameter which lead, after more than 
50 years of struggling, to a microscopic understanding of 
the phenomenon of superconductivity. 

BCS theory, however, cannot be considered a predic- 
tive theory in the sense that it would allow the compu- 
tation of material-specific properties. Moreover, mate- 
rials with strong electron-phonon coupling, such as nio- 
bium or lead, are poorly described by BCS theory. In 
these strong-coupling materials, phonon retardation ef- 
fects play a very important role. A proper treatment of 
those effects was developed by Eliashberg 2 ' 3 . His theory 
can be viewed as a GW approximation^ in terms of the 
Nambu-Gorkovai Green's functions. Eliashberg's theory 
not only achieves a successful description of the strong 
coupling simple metals like Nb and Pb, it also provides 
a convincing explanation of the superconducting features 
of more complex materials such as MgB^. 



In spite of its tremendous success, Eliashberg theory, in 
its practical implementation, has to be considered a semi- 
phenomenological theory. While the electron-phonon in- 
teraction is perfectly accounted for, correlation effects 
due to the electron-electron Coulomb repulsion are diffi- 
cult to handle in this theory. Those effects are condensed 
in a single parameter, /i*, which represents a measure 
of the effective electronic repulsion. Although /i* could, 
in principle, be calculated by diagrammatic techniques 3 , 
first principles estimates of /i* are extremely hard to 
make, and in practice, \i* is treated as an adjustable pa- 
rameter, usually chosen such that the experimental T c is 
reproduced. 

The goal of this work is to develop a true ab-initio 
theory for superconductivity which does not contain any 
adjustable parameters. The crucial point is to treat the 
electron-phonon interaction and the Coulombic electron- 
electron repulsion on the same footing. This is achieved 
within a density-functional framework. 

Density functional theory (DFT)L£*2i enjoys enormous 
popularity as an electronic-structure method in solid- 
state physics, quantum chemistry and materials science. 
DFT combines good accuracy with moderate numerical 
effort and is often the method of choice especially for 
large molecules and solids with a big unit cell. DFT is 
based on the Hohenbcrg-Kohni theorem which ensures 
a rigorous 1-1 correspondence between the ground-state 
density and the external potential. At finite temperature, 
the correspondence holds 10 between the density in ther- 
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mal equilibrium and the external potential. As a conse- 
quence, all physical observables of an interacting electron 
system become functionals of the density. The practical 
implementation of DFT rests on the Kohn-Sham& scheme 
which maps the interacting system of interest on an aux- 
iliary non-interacting system with the same ground-state 
thermal density. 

Traditional DFT, by its very nature, inevitably 
involves the Born-Oppenheimer approximation: One 
is supposed to calculate the electronic ground- 
state/thermal density that corresponds 1-1 to the elec- 
trostatic potential of clamped nuclei. To overcome this 
limitation Kreibich and Gross (KG)ii recently presented 
a multicomponent DFT which treats both electrons and 
nuclei quantum mechanically on the same footing. The 
KG theory involves two "densities" : The electronic den- 
sity n(r) referring to a body-fixed coordinate frame, and 
the diagonal T(R) of the nuclear iV-body density matrix. 
The exchange-correlation functional appearing in the re- 
sulting Kohn-Sham equations depends on both "densi- 
ties" and contains, formally, all non-adiabatic couplings 
between the electrons and the nuclei. 

Hence, in principle, the KG framework should be able 
to describe (conventional) superconductivity. In practice, 
however, it is not advisable to attempt such a description. 
The situation is quite similar to the DFT treatment of 
magnetic effects: By virtue of the Hohenberg-Kohn theo- 
rem, magnetic effects can be described on the basis of the 
density alone. In particular, the order parameter of spin 
magnetism, the spin magnetization m(r), is a functional 
of the density m = m[n]. However, this functional has 
to be highly non-local and its explicit form is unknown. 
Therefore, it is advisable to include the order parameter 
m(r) as an additional "density" in the DFT formula- 
tion— This version of DFT is known as spin-DFT. It 
is the standard form of DFT which is employed in all 
practical applications. 

A similar idea was suggested in 1988 by Oliveira, Gross 
and Kohn (OGK)H to treat superconductors in a DFT 
framework. OGK proposed the inclusion of the order pa- 
rameter x( r ) r ') that characterizes the superconducting 
phase as a basic "density" in the DFT formulation. OGK 
dealt with singlet order parameters. A generalization to 
triplet order parameters was given later—. 

The complexities of the many-body problem were cast 
into an exchange and correlation term, but in contrast 
to ordinary density functional theory where a variety 
of functionals has appeared over the past thirty years, 
very few exchange-correlation functionals have been pro- 
posed for the superconducting state. To our knowledge, 
only a local density approximation describing the purely 
electronic interactions has been presented— However, 
the usefulness of the OGK approach was demonstrated 
by Gyorffy and coworkers in their study of niobium and 
YBa2Cu307 using a semi-phenomenological parametriza- 
tion of the exchange-correlation functionaliL— 

The OGK formulation was triggered by the discovery 
of high-T c superconductors^ where an electronic pair- 



ing mechanism was believed to be dominant. Hence, the 
OGK description treats the Coulomb repulsion between 
electrons formally exactly while the electron phonon cou- 
pling only enters through a given, non-retarded BCS-type 
electron-electron interaction, i.e. strong electron-phonon 
coupling cannot be dealt with. 

In this work we develop a DFT for superconductors 
based on the three densities n(r), x(r, r') and T(R). The 
formalism can thus be viewed as a superconducting gen- 
eralization of KG or as a strong coupling generalization 
of OGK. It leads to a set of - formally exact - Kohn-Sham 
equations for electrons and nuclei. These equations con- 
tain exchange-correlation potentials which are universal 
functionals of the three densities n, x an d T. For the time 
being, we do not study effects found in the presence of 
magnetic fields. Those can be accomodated by general- 
izing the framework to include the current density as an 
additional variabl e 20 ' 21 . 

The success of any density functional theory crucially 
depends on the availability of accurate approximations 
for the exchange-correlation functionals. The main body 
of this article is devoted to the construction of such ap- 
proximate exchange-correlation functionals. Diagram- 
matic many-body perturbation theory is used for this 
purpose. In a second article (henceforth referred to as 
II), these approximate functionals are employed to cal- 
culate superconducting properties of elemental metals. 

The present paper is organized as follows: In Sect. [H] 
we derive a multicomponent DFT for the superconduct- 
ing state. This theory leads to a set of Kohn-Sham equa- 
tions that are described in the following section. Sect. II VI 
is devoted to the development of Kohn-Sham perturba- 
tion theory. The resulting exchange-correlation poten- 
tials are discussed in Sect. [V] A simple BCS-like model 
is described in Sect. I VII This model is used, in Sect. lVIll 
to analyze the approximate exchange-correlation kernels 
entering the linearized DFT gap equation. Finally, in 
Sect. IVIIII the exchange-correlation contributions to the 
non-linear gap equation are discussed. 



II. MULTICOMPONENT DFT FOR 
SUPERCONDUCTORS 

It is clear that a balanced treatment of the electron- 
phonon and Coulomb interactions has to start from the 
many-body electron-nuclear Hamiltonian (atomic units 
are used throughout this paper) 

H = f c + U cc + f n + U nn + U CD , (1) 

where T° represents the electronic kinetic energy, U ee the 
electron-electron interaction, T n the nuclear kinetic en- 
ergy and U nn the Coulomb repulsion between the nuclei. 
The interaction between the electrons and the nuclei is 
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described by the term: 



U c 



E 



d 3 R$l(r)&(R)-. 



R 



(2) 

where ^(r) and $>>(R) are respectively electron and 
nuclear creation operators. (For simplicity we assume 
the nuclei to be identical, and we neglect the nuclear spin 
degrees of freedom. The extension of this framework to 
a more general case is straightforward.) Note that there 
is no external potential in the Hamiltonian. 

To develop a multicomponent DFT for the electron- 
nuclear system we have to proceed with care. The 
Hamiltonian JQl describes a translationally invariant and 
isotropic system. Thus, both the electronic and nuclear 
one-particle densities are constant and therefore not use- 
ful to characterize the system. This problem can be 
solved by adopting a body-fixed reference fram e 11 ! 22 . In 
this article we are interested in infinite solids where the 
nuclei perform small oscillations around the equilibrium 
positions. Furthermore, we assume that the solid is not 
rotating as a whole. Fortunately, in this case, the body- 
fixed reference frame coincides with the normal Cartesian 
system commonly used to describe solids. The situation 
is very different for finite systems, which have to be han- 
dled by using appropriate internal coordinates. 

In order to formulate a Hohcnberg-Kohn type state- 
ment, the Hamiltonian Q is generalized to 

H = f c + f a + U cn + U cc + V c c xt + V c a xt + A oxt -fiN. (3) 

The external potential for the electrons is defined as 



Kxt = E 



d 3 r¥ <T (r)v e ext (r)* a (r). 



(4) 



Since, at this level, the nuclei are taken into account ex- 
plicitly, the lattice potential is not treated as an external 
field, but is included via the interaction term U cn . The 
term V* xt is introduced as a mathematical trick to prove 
the Hohenberg-Kohn theorem, and will be taken to zero 
at the end of the derivation, t^xt i s a multiplicative N- 
body operator with respect to the nuclear coordinates 



(5) 
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V£t = Jd 3 Rv^(R)T(R), 

where we have defined R = {Ri, R2, • • • , Rn}, d 3 R 
d 3 R 1 d 3 R 2 ---d 3 R N , and 

f (B) = *t( fll ) . . . &(R N )$(R N ) . . . ( 6 ) 

is the diagonal part of the nuclear TV-particle density ma- 
trix operator. Note that the term V^t includes the inter- 
action between the nuclei U nn (to which it reduces if no 
other external nuclear potentials are present). The term 



A nx t — 



dV 



describes an external pairing field, and usually vanishes 
unless our system is in the proximity of an adjacent su- 
perconductor. However, this term is required to break 
the gauge invariance of the Hamiltonian, and the limit 
A C xt - ► can only be taken at the end of the deriva- 
tion. Note that the term 10 describes a singlet pair- 
ing field. The extension to triplet superconductors is 
straightforward 14 ;. Finally, /x stands for the chemical po- 
tential, and N is the number operator for the electrons 
(we treat the electronic degrees of freedom in a grand- 
canonical ensemble). 

Our multicomponent formulation is based on three 
densities: 



i) The electronic density 



(8) 



which is defined in the usual way. The bracket (• • • ) 
denotes the thermal average (A) = Tr{poA}, with 
the grand canonical statistical density operator p a — 
c^^/Trjc - ^ } in the superconducting state. We fur- 
thermore define the inverse temperature (3 = 1/T. 

ii) The anomalous density 



X{r,r') = O&tM*; 



(9) 



(7) 



which is the order parameter characterizing the singlet 
superconducting state. This quantity is finite for super- 
conductors below the transition temperature and zero 
above this temperature. 

iii) To describe the nuclear degrees of freedom, we use 
the diagonal part of the nuclear iV-particle density matrix 
T(R) = (T(R))- Alternatively, one could define a mul- 
ticomponent DFT using the one-particle density for the 
nuclei n n (R) = (<t'(R)<&(R)). However, in the following 
it will be convenient to transform the nuclear degrees of 
freedom to collective (phonon) coordinates. Using n n (R) 
would lead to a one body equation for non-interacting nu- 
clei. Thus, the nuclear Kohn-Sham equation would not 
lead to realistic phonons with a proper dispersion rela- 
tion. Only Einstein phonons could be present in this 
system. This is also clear from the fact that a system 
of non-interacting particles does not exhibit collective 
modes. With our choice of T(R), the nuclei obey a TV- 
particle equation, very similar to the Born-Oppcnhcimcr 
equation, and where phonon coordinates can be easily 
introduced. 

As usual, the Hohenberg-Kohn theorem guaran- 
tees a one-to-one mapping between the set of 
the densities {n(r), x(r,r'),r(R)} in thermal equi- 
librium and the set of their conjugate potentials 
{vl xt (r), A ext (r, r'), v^ xt (R)}. As a consequence all ob- 
servables are functionals of the set of densities. Finally 
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it assures that the grand canonical potential, 

n[n,x,r]=F[n, X ,T] + J d 3 rn(r)[v° xt (r) — fj] 
r d» r /dV W ,,OA^(,,/) + /,.c] 

+ |d 3 Rr(RK' x ,(B), (io) 

is minimized by the equilibrium densities. We use the 
notation A[f] to denote that A is a functional of /. The 
functional F[n, Xi T] is universal, in the sense that it does 
not depend on the external potentials, and is defined by 

F[n, X ,T] =T>, X ,r]+T>, X ,r] 

+ c/°>, x , r] + u cc [n, x , r] - ^s[n, x, r] (11) 

where 5 stands for the entropy of the system 

S[n, X , T] = -Tv{p [n, x, r] log(p [n, x, r] )} . (12) 

The proof of the theorem follows closely the proof of the 
Hohenberg-Kohn theorem at finite temperatures^ and 
will not be presented here. 

In standard DFT one normally defines a Kohn Sham 
system, i.e. a non-interacting system chosen such that 
it has the same ground-state density as the interacting 
one. In our formulation, the Kohn-Sham system consists 
of non-interacting (superconducting) electrons, and in- 
teracting nuclei. It is described by the thermodynamic 
potential [cf. Eq. JTJJ] 

n s [n,x,r] = F s [n, x ,T] + J d 3 r n(r)[<(r) - /i 8 ] 

+ J d"R r(R)vf(R) , (13) 

where F s if the counterpart of 111(1 for the Kohn-Sham 
system, i.e., 

F B [n, X ,T] =T>,x,r]+T s >,x,r]-is s [n,x,r]. (14) 

P 

Here T°{n, x, T], T s n [n, x, T] and S s [n,x,r] are the elec- 
tronic and nuclear kinetic energies and the entropy of 
the Kohn-Sham system, respectively. 

From Eq. ((13(1 it is clear that the Kohn-Sham nuclei 
interact with each other through the iV-body potential 
vf(R) while they do not interact with the electrons. 

By applying the Hohenberg-Kohn theorem to both the 
interacting and the non- interacting systems, and requir- 
ing the densities of the Kohn-Sham system to reproduce 
the densities of the fully interacting one, we can identify 
the expressions for the effective Kohn-Sham potentials. 
As usual, these include contributions from external fields, 



Hartree, and exchange-correlation terms. The latter ac- 
count for all the many-body effects stemming from the 
electron-electron and electron- nuclear interactions. To 
simplify the expressions, we now set the auxiliary exter- 
nal potentials to zero. 

The Kohn-Sham potential for the electrons v°(r) reads 

as 



v s [n, x 



^{r) = -ZY^J^R- 



nm 



\r — r 



-+<>,x,r](r). (15) 



The first term, the electron-nuclear Hartree potential, 
reduces to the usual nuclear attraction potential if we 
assume that the nuclei are classical and perfectly local- 
ized at their equilibrium positions. This term is usu- 
ally treated as the external potential in standard DFT. 
The last two contributions to v*(r) are respectively the 
Hartree potential, which accounts for the classical repul- 
sion between the electrons, and the exchange-correlation 
term. 

The anomalous Kohn-Sham potential A s is given by 
A s [n, x, r] (r, r') = + AxcK X , r] (r, r') . (16) 

Note that the first term, the so-called anomalous Hartree 
potential, gives rise to a positive contribution to the en- 
ergy. 

Finally, the nuclear potential is 



Vs >,x,r](£) = ]T 



z 2 



Ra — R/3 



n(r) 



r , - p | +<>,X,r](H). (17) 
r - R a \ 



The first term stems from U nn , and describes the bare 
nuclear-nuclear repulsion. The second is the nuclear- 
electron Hartree term and is the counterpart of the first 
term in Eq. ((15(1 . 

As in standard DFT, the exchange-correlation poten- 
tials are defined as functional derivatives 

< c [n,x,r](r) = ^^' r] , (18a) 
A xc [n,x,r](r,r') = - ^M , (18b) 



v» c [n,x,T](E) 



5x*(r, r') 
5F xc [n, X ,r] 
ST(R) ' 



(18c) 



The exchange-correlation free energy is defined through 
the equation 

F[n, X ,T] = F s [n,x,r]+F xc [n,x,r] 

+ U^[T}+E^[n,x}+E^[n 7 r}. (19) 
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There are two contributions to E^, one stemming from 
the electronic Hartree potential, and the other from the 
anomalous Hartree potential 



*%M = \ fd 3 rjd 3 r 



, n(r)n(r') 
\r — r'\ 

d'r/dV |a f (r>r T - (20) 



Finally, denotes the electron-nuclear Hartree energy 



n(r)T(R) 



R a 



(21) 



III. THE KOHN-SHAM EQUATIONS 

The problem of minimizing the Kohn-Sham grand 
canonical potential (|13|l can be transformed into a set 
of three differential equations that have to be solved self- 
consistently: Two of them are coupled and describe the 
electronic degrees of freedom. Their algebraic structure 
is similar to the Bogoliubov-de Gennes equations. The 
third is an equation for the nuclei resembling the familiar 
nuclear Born-Oppenheimer equation. 



A. Electronic equations 

The Kohn-Sham Bogoliubov-de Gennes (KS-BdG) 
equations*^ read 



+ »»(»•) - A* 



V 2 

-— +vl(r)-n 



Ui(r) + J dV A s (r,r')vi(r') 

= EiUi(r), (22a) 

Vi(r) + /dV A B *(r,r')ui(r') 

= E iVi (r), (22b) 



where Ui(r) and Vi(r) are the particle and hole ampli- 
tudes. This equation is very similar to the Kohn-Sham 
equations in the OGK formalism 13 . However, in our for- 
mulation the lattice potential is not considered as an ex- 
ternal potential but enters via the electron-ion Hartree 
term. Furthermore, our exchange-correlation potentials 
depend parametrically on the nuclear density matrix, and 
therefore on the phonons. 

Although these equations have the structure of static 
equations, they contain, in principle, all retardation ef- 
fects through the exchange-correlation potentials. 

A direct solution of the Kohn-Sham Bogoliubov-de 
Gennes equations-^ is faced with the problem that one 
needs extremely high accuracy to resolve the supercon- 
ducting energy gap, which is about three orders of mag- 
nitude smaller than typical electronic energies. At the 
same time, one has to cover the whole energy range of 



the electronic band structure. The so-called decoupling 
approximation 2 ^^ 2 ^ relieves the problem by separating 
these two energy scales. 

The particle and hole amplitudes can be expanded in 
the complete set of wave functions {y>i} of the normal- 
state Kohn-Sham equation 



(fi(r) = ei cpi(r) (23) 



which can be solved by standard band-structure meth- 
ods. The decoupling approximation then implies the fol- 
lowing form for the particle and hole amplitudes, where 
only one term of the expansion is retained: 



Ui(r) « Ui(fi{r) ; Vi(r) w vnpi(r) . 



(24) 



In this way the eigenvalues in Eq. 122H become E t = ±Ei, 
where 



i^^ + IA.I 2 , (25) 

and £i — ti — [i. This form of the eigenenergies allows 
the interpretation of the matrix elements of the pair po- 
tential as the gap function of the superconductor. The 
coefficients Ui and Vi also have simple expressions within 
this approximation 



U, 




Finally, the matrix elements Aj are defined as 
A, : /d 3 rfd 3 rV:(r)A s (r,r')^(r'), 



(26a) 
(26b) 

(27) 



and fa is the phase e 1 *^ = Aj/|Aj|. Within the decou- 
pling approximation, the normal and the anomalous den- 
sities can be easily obtained from the particle and hole 
amplitudes 



i 

X{r,r') = 1^2 



1- Jitanhf^ 
Ei I 2 



|^(r)| 2 (28a) 



^ta n h(|s < ) r -,(r)^(r'. ) .2Xh, 



The validity and limitations of the decoupling approx- 
imation will be discussed in detail in II. 



B. Nuclear equation 

The Kohn-Sham equation for the nuclei has the form 

*iCH) = 5*iGR)- (29) 



2M 
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This equation has the same structure as the usual nu- 
clear Born-Oppcnheimer equation. We emphasize that 
the Kohn-Sham equation 1)29(1 does not rely on any ap- 
proximation and is, in principle, exact. As already men- 
tioned, we are interested in solids at relatively low tem- 
perature, where the nuclei perform small amplitude oscil- 
lations around their equilibrium positions. In this case, 
we can expand the v"[n, x,T] in a Taylor series around 
the equilibrium positions, and transform the nuclear de- 
grees of freedom into collective (phonon) coordinates. 
In harmonic order, the nuclear Kohn-Sham Hamiltonian 
then reads 



\.q 



(30) 



where &x,q are the phonon eigenfrequencies, and b\^ q de- 
stroys a phonon from the branch A and wave-vector q. 
Note that the phonon eigenfrequencies are functionals of 
the set of densities {n,x,r}, and can therefore be af- 
fected by the superconducting order parameter. Within 
the harmonic approximation, the nuclear density matrix 
T(R) is given by 



T(R) = J2MVx, q )\hx, q (Q)\ 2 

X,q 



(31) 



where np{VL) denote the Bose occupation numbers and 
hx,q{Q) arc harmonic oscillator wave functions referring 
to the collective coordinates Q. 



C. Gap equation 

In Fig. ^ we sketch the Kohn-Sham self-consistent pro- 
cedure within the decoupling approximation. We start by 
finding suitable approximations for the Kohn-Sham po- 
tentials to start the cycle: for w°'°[n,x,r] we can use the 
Kohn-Sham potential stemming from a standard DFT 
calculation for the non-superconducting ground-state, i.e, 
v^ s [n] . This is a very good approximation as the depen- 
dence of f|[n, x,r] on x an d T is certainly very weak 
for the usual superconductors at low temperature; The 
starting pair potential Ag[n,x,r] can be approximated 
by a square well centered at the Fermi energy with width 
of the order of the Debye frequency and height computed 
from a BCS model; Finally, for u" ,0 [n,x,r] we can use 
the ground-state Born-Oppenheimer potential. The next 
two steps of the self-consistent cycle can be performed in 
parallel: 

• Equation (|23|l is solved to obtain the wave- 
functions c/?;'s and the eigenenergies e^'s. These can 
then be used within the decoupling approximation, 
equation (|24[1 . to calculate the normal and anoma- 
lous densities through Eqs. I|28|) . We note that the 
chemical potential fi entering Eq. (|23|l has to be 
adjusted at every iteration, such that the density 
n(r) integrates to the correct particle number N. 



Approximate 
<>, X ,r](r) , A°[n, x ,r](r,r') , X , T}(R) 



Solve electronic KS 
Eq. (23) for tpi(r) 



Solve nuclear KS 
Eq. (29) in the 
harmonic approximation 



decoupling approximation 



calculate n, x 
through Eq. (28) 



calculate T 
through Eq. (31) 







calculate th 

v%[n,x 
A s [n,x, 
<[n,x 


e potentials 

,r](r), 

r](r,r'), 

,r]GR) 



self-consistency 
convergence check 



No 



output 



FIG. 1: Schematic flow chart for the iterative Kohn-Sham 
scheme within the decoupling approximation. 



• With w"[n, x, r] we solve the nuclear equation 129|) 
by expanding the nuclear potential to harmonic 
order to obtain the phonon eigenfrequencies and 
eigenmodes. The nuclear density matrix T then 
follows from Eq. (j^l . 

Finally, the set of densities {n, x, T} is used to evaluate 
the new Kohn-Sham potentials {vf , A s , w"} from the def- 
initions I|15ll6ll7fl . At this point, if self-consistency is 
reached, the cycle is stopped. Otherwise, the new poten- 
tials are used to restart the cycle. 

It is clear that, even within the decoupling approxi- 
mation, the task of solving the self-consistent cycle de- 
picted in Fig. ^ is rather demanding. Furthermore, we 
are required to provide (good) approximations for the 
three functionals vf[n,x,r], A s [n, x, T], and w"[n, x,T]. 
As our objective is to study superconductivity, we will 
make two simplifying assumptions: i) Ug[n, x,T] can be 
well approximated by the ground state functional used 
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in standard density functional theory Ug [n]. As the 
energy scale of the phonons and of the superconduct- 
ing gap is three orders of magnitude smaller than elec- 
tronic energy scales (like the Fermi energy) this is a 
very reasonable assumption, ii) The nuclear functional 
Vg[n, Xi r] is approximated by the ground-state Born- 
Oppenheimer energy surface. It is well known that calcu- 
lations performed within the Born-Oppenheimer approx- 
imation yield phonon frequencies that are in excellent 
agreement with experimental dispersions^ 6 -. We therefore 
expect this to be an excellent approximation to the Kohn- 
Sham phonons. However, we are neglecting the effect of 
the superconducting pair potential on the phonon dis- 
persion. This effect has been measured experimentally^, 
and it turns out to be quite small. Note that these two 
approximations are equivalent to fixing Vg[n, and 
vf[n, X ,T] to u s e '°[n, X ,r] and r]. 

By inserting Eqs. I|28|) in Eq. (|16[) and by subsequently 
inserting Eq. (|16[) on the right-hand side of Eq. I|27l) , we 
obtain the DFT gap equation 



A,; = A 



Hxc i 



[Mi Aj 



(32) 



where Ah xc stands for the sum of the Hartree and 
exchange-correlation contributions to the functional. 
Note that through the exchange-correlation functional 
the right-hand side of Eq. IMl'l) becomes a highly com- 
plicated functional of fi and A^ . The dependence on the 
gap function is totally different from the usual mean-field 



gap equation (cf. Sect IVIIIjl . 

After these simplifying approximations, a Kohn-Sham 
calculation proceeds as follows: i) Using a standard 
band structure code, the ground state wave- functions and 
band structure are obtained, ii) The Born-Oppenheimer 
phonon frequencies and eigenmodes are obtained from 
linear-responseSi calculations, again using standard tools 
widely available to the community, iii) The gap equa- 
tion lj3~2l is iterated until self-consistency is achieved. We 
can now see how the different energy scales are sepa- 
rated: The normal density, the anomalous density, and 
the phonon properties are obtained from three separate 
equations. 

In the vicinity of the transition temperature, the 
anomalous density will be vanishingly small. In this 
regime, the gap-equation can be linearized in x, leading 
to 



1 



^ -^Hxc i.j \p] 



l au.1.1 ( 



A, 



(33) 



where the anomalous Hartree exchange-correlation kernel 
of the homogeneous integral equation reads 



XC l.J 



[/'] 



SA 



Hxc i 



S X^ S Xj 



(34) 

We emphasize that the linearized gap equation can only 
be used to calculate the superconducting transition tem- 
perature. In particular, the function A^ that stems from 



the solution of this equation does not have any physical 
interpretation. 

We can gain further insight into Eq. I|33|l if we separate 
the kernel JFhxc i j into a purely diagonal term, Zi, and 
a non-diagonal part, K-i.j 



Aj = -Zi\ix] Aj - K, itj \p\- 



tanh 



A, , (35) 



This equation has the same structure as the BCS gap 
equation with the kernel K% j replacing the model interac- 
tion of BCS theory. On the other hand, Zi plays a similar 
role as the renormalization term in the Eliashberg equa- 
tions. This similarity allows one to interpret the kernel 
JCij as an effective interaction responsible for the binding 
of the Cooper pairs. The function ICij is a quantity of 
central importance in the density functional theory for 
superconductors. By studying fC^j for a specific mate- 
rial as a function of i and j one can learn which orbitals 
are responsible for superconductivity and, ultimately, by 
identifying those parts of the exchange-correlation func- 
tional (phononic/Coulombic) that lead to an effective at- 
traction, one can trace the mechanism responsible for the 
superconducting state. 

Note that Eq. (|35|l is considerably simpler than the 
Eliashberg equations as it does not dependent explicitly 
on the frequency. However, phonon retardation effects 
are included through the exchange-correlation terms. 
Furthermore, Eq. ll-iol) is not a mean-field equation like in 
BCS theory but contains correlation effects. A linearized 
gap equation can also be derived without the decoupling 
approximation^ 2 ^, leading to a similar equation, but 
with a four-point kernel. From this point of view, the 
decoupling approximation can be viewed as a diagonal 
approximation to this four-point kernel, neglecting the 
corresponding off-diagonal elements. 



IV. KOHN-SHAM PERTURBATION THEORY 

In the previous sections we showed how to develop a 
density functional theory for the superconducting state. 
The main equation of this theory, the gap equation i|32[l , 
allows, in principle, the calculation of the superconduct- 
ing gap for any system. However, to solve this equa- 
tion one needs approximations for A xc , the exchange- 
correlation contribution to the Kohn-Sham pair poten- 
tial. In the following, we will develop such approxima- 
tions by applying Kohn-Sham perturbation theory, as de- 
scribed by Gorling and LevyS, to superconducting sys- 
tems2^21i. This perturbation theory, that will treat both 
the electron-electron and electron-phonon interactions on 
the same footing, is a generalization of the method used 
by Kurth et al. to calculate the exchange-correlation en- 
ergy of the uniform superconducting electron gas^. 

Our starting point is the Hamiltonian of the electron- 
nuclear system as defined in Eq. @. This Hamiltonian 
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is then split into a suitably chosen reference Hamilto- 
nian H and the remainder, which is treated as a per- 
turbation. The most appropriate reference system for 
this formalism is as follows: i) The nuclear Kohn-Sham 
Hamiltonian (|29() rigorously defines the nuclear equilib- 
rium positions Rq. When expanded to harmonic order 
around these positions it can be written as the phonon 
Hamiltonian H ph (|5Ujl. ii) Next we define a Born- 
Oppenheimer Kohn-Sham Hamiltonian via a rigid-lattice 
potential given by the equilibrium coordinates R , 



H 



BO 



T c + V 



latt ,Rq 



Vr 



Hxc 



A 



Hxc 



(36) 



where 



d d r 



, n(r>) 



— + v c (r) 
\r-r'\ ^ Uxc[ '> 



(37) 



while Ahxc includes the anomalous Hartree and 
exchange-correlation contributions 



A H xc = - /d 3 r/dV * T (r)^ x (r') 



X*(r,r') 



H.c. 



(38) 



With the choice Hq = H ph + Hqq, the interaction 
Hamiltonian reads: 

H, = L> ec + L>° G ph - V£ xc - V£ xc - Ahxc (39) 

The Born-Oppenheimer electron-phonon coupling oper- 
ator £7gQ Ph is given by 



e— ph 



BO 



EE /d 3 rF A B ? °(r)*t(r)i(r)$ A , ? . 



(40) 



where V^°(r) is basically the gradient of the electronic 
Kohn-Sham potential with respect to the nuclear co- 
ordinates and the phononic field operator is <3E>a,<j = 
Note that now we have set the auxiliary 



3A, q 



'-\,- r 



external potentials V° Kt and A cxt to zero, and V^ t to the 
bare intcrnuclcar interaction. 

We believe that this is the most physical way to split 
the Hamiltonian, since the electronic structure calcula- 
tion for n(r), in practice, is usually performed for fixed 
nuclear positions; the nuclear dynamics is absorbed in 
the exchange correlation functionals. Furthermore, stan- 
dard calculations for the electron-phonon coupling, which 
are based on linear response theory, involve the above 
coupling potentials V^°(r). Note that, besides the in- 
teraction terms C/ ee and U^q 1 , the perturbation includes 
the Hartree and exchange-correlation potentials. In ap- 
pendix [S] we give some more details of this construction. 



We now develop a many-body perturbational ap- 
proach, which will ultimately lead to explicit expressions 
for the exchange-correlation functionals. The construc- 
tion of this approach follows closely the standard many- 
body perturbation theory described in many textbooks^. 
Our objective is to expand the difference Ail = £1 — fi s 
in a power series. From this difference and with the def- 
initions I|1U|) and l|13|) it is straightforward to derive an 
expression for the functional F xc . 

In standard perturbation theory Af2 is written as an 
expansion in powers of e 2 and g 2 , where e (the electron 
charge) and g (the electron-phonon coupling constant) 
measure respectively the strength of the Coulomb and of 
the electron-phonon interactions. In our approach, how- 
ever, every order in perturbation theory contains all or- 
ders in e 2 and g 2 . This is due to the special form of the 
perturbation Hi that involves the exchange-correlation 
potentials which contain implicitly all orders in the inter- 
actions. Therefore, for book-keeping purposes, we mul- 
tiply the perturbation H± by a dimcnsionless parameter 
A. In this way, the terms appearing in the perturbative 
expansion can be labeled by powers of A. 

The grand canonical potential f2 can be written as 



n = ~iog(z), 



(41) 



where the partition function has the definition Z = 
Tr{exp(— I3H)}. From this expression it follows directly 
that 



fi-fi s = -Ilog(| 



(42) 



where Z s is the partition function of the Kohn-Sham sys- 
tem. It is then straightforward, using the standard ma- 
chinery of many-body perturbation theory, to write the 
ratio Z /Z s as a series expansion in A which can be eval- 
uated with the help of Wick's theorem. Moreover, the 
number of terms in the series can be reduced by using a 
generalization of the linked-cluster theorem^. The final 
result reads 



fi — il s = — — E]( a U connected diagrams) . 

P 



(43) 



In this expression the sum runs over all connected Feyn- 
man diagrams that are topologically distinct. (Alterna- 
tively, one can expand diagrammatically the propagators 
and use the Galitskii-Migdal formula to evaluate the en- 
ergy^. The two approaches are equivalent.) 

There are several Kohn-Sham propagators that enter 
the Feynman diagrams: First we have the contraction 
that reduces to the usual Green's function for systems 
that are not superconducting 



G s CTCT ,(rr,rV) = -<f ^(rr)^, (rV)>, , 



(44) 



where T is the time-ordering operator, which orders the 
operators from right to left in ascending imaginary time 



order, and the average (• • • ) s is done with respect to the 
Kohn-Sham statistical density operator p s . This Green's 
function is represented in the Feynman diagrams by a line 
with two arrows pointing in the same direction. Further- 
more, due to the presence of pairing fields in the Kohn- 
Sham system - the following (anomalous) averages 
are non- vanishing for superconducting systems 




FIG. 2: Lowest order electronic (a) and phononic (b, c) con- 
tributions to F xc . 



F^,(rr,rV) = -<f ^(rr)Vv (rV)) 8 , (45a) 
F s J,(rr,r'r') = -<f & (rr)#,(rY)> B . (45b) 



In Feynman diagrams these propagators appear as lines 
with two arrows pointing outwards (for F s ) and as 
lines with two arrows pointing inwards (for F s t). (The 
Green's function (|44|l and the anomalous averages H45|) 
can be conveniently assembled into a matrix Green's 
function in Nambu-Gorkov spaced.) Finally, as the per- 
turbation includes the electron-phonon interaction term 
-ffc-ph, some diagrams contain the phonon propagator 
(represented as a curly line) 

Dl q (r,r') = {U Kq (r)^l q (r')) a . (46) 

As usual in finite temperature many-body theory, it is 
convenient to work in imaginary frequency space. The 
time Fourier transform of the Green's function (|45|l is 
defined as 

G^,(rr,rV) ^Se-^GJ^r.r'jw,,), (47) 

where u n = (2n + 1)tt / j3 are the odd Matsubara frequen- 
cies. The frequency dependent anomalous propagators 
have a similar definition. In Matsubara space the Kohn- 
Sham propagators read, in terms of the Kohn-Sham par- 
ticle and hole amplitudes and of the Kohn-Sham eigenen- 
ergies, 



E 



Ui (r)u*(r') Vi(r)v*(r') 



E,. 



10,',, 



(48a) 



F*y(r,r';cj n ) = 5„ t s8gp.(</) 

~v*{r)u t {r') Ul {r)v*{r') 



\u> n — Ei 



(48b) 



F^,(r,r';cj„) = d CT _ CT ,sgn(o-) 

~u*(r)vi(r') Vl {r)u*{r') 



E 



iuj n + Ei 



lUn — Ei 



(48c) 



On the other hand, the phonon propagator depends on 
the even Matsubara frequencies v n = 2n7r//3, 



DIM 



,,2 I Q 2 ' 



(49) 



In first order in A there is only one diagram contribut- 
ing to F xc . This diagram, depicted in Fig.|2i, is of purely 
electronic origin and has the same form as the standard 
exchange diagram. (The wavy line in Fig. [2b represents 
the Coulomb interaction.) This term can be written as 
(for simplicity we write all contributions to F KC within 
the decoupling approximation 



1 I tanh U Ei 



1 - 



E, 



tanh | 



, (50) 



where the matrix elements of the Coulomb interaction 
are defined as 



J %,3 



= d 3 r dV^(r) W (r') 



1 



T <p {r)ip*{r') . 



(51) 

As the expectation value (&\, q ) — 0, the electron-phonon 
interaction does not contribute to -F xc in first order per- 
turbation theory in A. The first non- vanishing terms 
appear in second order in A (first order in g 2 ) and are 
depicted in Fig. [2b, c. The first of these terms (Fig. [2b) is 
the counterpart of the anomalous term that contributes 
to the electronic Hartree energy (|2(J[) . Its analytic form 
can be written as 



A,q ij 



2 A, A* 



EiEj 



I(Ei,Ej, SlA,q) 



-I(E t ,-E 3 ,n x , q )\, (52) 

where we defined the matrix elements of the electron- 
phonon coupling constant 



9{ q = /d^:(r)F A B >V,(r), (53) 



while the function / is 

I(E,E,Sl) = J2 ■ -Eiuj 2 - E' Jul - oj 2 ) 2 + « 2 ' 

(54) 

The three fractions in the sum come from the denomina- 
tors of the two Green's functions G s and from the phonon 
propagator D B . It is possible to perform the frequency 
sums using standard complex contour integration meth- 
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ods. The final result is 



I(E i ,E j ,Sl) = f p (E i )f l3 (E J )np(n) 



Ei — Ej 



Ei — Ei 



n 



(55) 



A diagram analogous to the one depicted in Fig. [2b 
but with the phonon propagator replaced by the bare 
Coulomb interaction exists as well. This diagram is the 
anomalous Hartree term which appears as the second 
term on the right hand side of Eq. (|2U[1 . Since this term is 
treated explicitly in the electronic Kohn-Sham equations 
it is not part of the exchange-correlation functional. 

Note that expression l|52l) is so much more compli- 
cated than the anomalous Hartree term simply because 
the phonon propagator describes a retarded interaction. 
If we assumed a retarded electronic interaction instead 
of the instantaneous Coulomb potential l/\r — r'\ the 
anomalous Hartree contribution would look very similar 

to gg&. 

The second phononic term that contributes to F xc is 
depicted in Fig. [3:. This Feynman diagram has the same 
structure as the electronic exchange term (Fig. Efe-) • How- 
ever, and again due to retardation effects, its analytic 
structure is more complicated than (|50|) . namely 



A,<J ij 



1+ J| In i[e - e > ■" a " ) 



(56) 



The self-energy diagrams contributing to Fs$ and F x ^ 
resemble the self-energy diagrams that appear in Eliash- 
berg theory 2 ^. By virtue of Migdal's theoremM, vertex 
corrections should be small and are therefore neglected, 
both in Eliashbcrg theory and in our treatment. There 
is, however, an important difference: in Eliashberg the- 
ory the propagators that enter the self-energy diagrams 
are dressed propagators, while in our case we have bare 
(Kohn-Sham) propagators. By using the bare propaga- 
tors we neglect more diagrams than those containing ver- 
tex corrections. We postpone a more detailed discussion 
of this problem to the section on the exchange-correlation 
kernels. 



FUNCTIONAL DERIVATIVES AND THE 
CHAIN RULE 



We have seen in Eq. (|18bl) how the anomalous 
exchange-correlation potential is defined as the func- 
tional derivative of the exchange-correlation free-energy 
functional with respect to the anomalous density %. How- 
ever, the contributions to the exchange-correlation free- 
energy functional that stem from Kohn-Sham pertur- 
bation theory are only known as explicit functionals of 



the Kohn-Sham orbitals (pi(r), the Kohn-Sham single- 
particle energies (e$ — /it), and the pair potential A*. Of 
course, the Hohenberg-Kohn theorem tells us that these 
quantities are themselves functionals of the densities, so 
the free-energy is an implicit functional of the densities. 
Furthermore, if one makes the additional approximation 
that i>x C does not depend on x then the Kohn-Sham or- 
bitals ifi(r) and the eigenenergies ti are also independent 
of the anomalous density. F xc is then a function of the 
chemical potential \i and a functional of the (complex) 
pair potential 



F xc = F xc [n, | A, 



(57) 



For convenience, instead of working with A^ and A*, 
we prefer using the modulus square of the pair potential 
and its phase as independent variables. The anomalous 
exchange-correlation potential can thus be evaluated us- 
ing the chain rule for functional derivatives 



A, 



SF XC 5^ 
Sfi Sx* 



■E 



SF XC 8\A 3 \ 2 SF XC Sfa) 



i) 5X* 



(58) 

The partial derivatives of F xc can now be calculated di- 
rectly. The remaining functional derivatives are somehow 
harder to obtain, but can be derived from the definitions 
of the densities, equations l|28Jl . and from the fact that 
the particle density and the anomalous density are inde- 
pendent variables. This latter condition can be expressed 
through the relation 



6n(x) 



0. 



Sx*(r,r') 

Moreover, we make use of the two trivial conditions 



(59) 



*x* 



Si, 



sx* 



0. 



(60) 



From the above conditions, and after some algebra, we 
arrive at the expressions for the remaining functional 
derivatives 



S\A< 



s_M 

SX* ' 



Yf 



iS 



E, 



Y 1 lA-l 2 — 

1 1 jl s x * 



" A* tanh (| ^) 



(61a) 
(61b) 

(61c) 



The functions Zf and Z} have the definitions 



Z? = 



Ei 2 tanh 



(fa) 



YP 



cosh' 



w 



(62) 
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and, finally, the two functions K° and read 



Y 



y.i = k. tanh ( —E, 

Ei 



^- tanh ( -Ei , 
^ V2 ; CO sh 2 (f^) 



2 « 



cosh 2 



(63a) 
(63b) 



approach is closer to many-body perturbation theory, 
and provides a natural relationship between the Green's 
function and the exchange-correlation potentials of DFT. 
However, both approaches lead to the same result if the 
free energy in the first method is expanded up to the same 
order in perturbation theory as the Green's functions in 
the second method. 



There exists another method to obtain exchange- 
correlation functionals using Kohn-Sham perturbation 
theory without resorting to functional derivatives. This 
method follows the ideas of Sham and Schliiten^, and 
provides a very direct connection between many-body 
perturbation theory and the normal and anomalous 
exchange-correlation functionals. In the following we will 
give a brief account of the Sham-Schliiter method for su- 
perconductors. 

There is a simple connection between the electron den- 
sity n(r) and the interacting many-body Green's function 



n(r)= limj^^e^Cfv;^). (64) 



The definition of the interacting Green's function is sim- 
ilar to Eq. H44[). but with the thermal average weighted 
by the interacting statistical operator po. The infinitesi- 
mal rj is used to ensure the correct ordering of the field 
operators in Eq. (fHfl so that Eq. (|54"Jl is fulfilled. In the 
same way it is simple to prove that the anomalous density 
obeys the relation 



X(r,r') 



lim — 



e"'"' 



F n (r,r';-w„), (65) 



where F is the interacting anomalous propagator. On 
the other hand, we defined the Kohn-Sham system as the 
non-interacting system whose both normal and anoma- 
lous densities are equal to the densities of the interacting 
system. Therefore, one can equally calculate the inter- 
acting densities from the Kohn-Sham Green's functions 



n(r) = lim — 

J?— 0+ f3 



EE ei " 



'Gw(r,r;w„) . 



(66) 



with an equation similar to Eq. I|65|l relating \ an d 
We then expand perturbatively the interacting Green's 
functions in l|64|) and Ij65(l . and equate the two equa- 
tions for n(r), and the two equations for x( r i r ')- As the 
perturbation (|39|l includes the terms V£ c and A xc , the re- 
sulting equations form a system of two coupled integral 
equations that allow the determination of i>° c and A xc . 
Those integral equations are the so-called Sham-Schliiter 
equations. 

The two methods to obtain the exchange-correlation 
potentials are conceptually quite different. The first uses 
the definitions l|18|l to derive the exchange-correlation po- 
tentials using a series of chain rules. The Sham-Schliiter 



VI. A SIMPLE BCS-LIKE MODEL 

We now introduce a simple model that will allow us 
to study in detail the functionals developed in this arti- 
cle. For simplicity, we assume that the pair potential 
has s-wave symmetry and behaves approximately like 
Ai = A(£j). This assumption is fulfilled by all tradi- 
tional superconductors. In this model, we can transform 
the gap equation into a one-dimensional equation in en- 
ergy space 

1 roo tanh ( f 

-2/ K'WMM) — ^-A(0, (67) 

where N(£) is the density of states at the energy £. It is 
possible to further simplify this equation by assuming a 
BCS-like model for both JC and Z. If we assume that the 
kernel JC and the renormalization term Z are constant 
in a shell of width lo c around the Fermi energy and zero 
outside this region, Eq. I|67|l can be solved analytically 
for the transition temperature T c 



T c ex exp 



1 + 2(0) 



AT(0)/C(0) 



(68) 



where JC(0) and Z(0) are the values of /C(£, £') and Z(£) 
at the Fermi surface. Equation (|68|l has exactly the same 
structure as McMillan's formula™^!, which is an approx- 
imate solution of the Eliashberg equations. This latter 
formula reads 



±r = — - exp 
1.20 v 



1.04(1 + A) 



A-/i*(l + 0.62A) 



(69) 



The number /it* , the Coulomb pseudo-potential of Eliash- 
berg theory, measures the strength of the electron- 
electron interaction. This parameter is quite hard to 
calculate and is often fitted to experimental data. As 
fi* is normally positive, it tends to decrease the super- 
conducting transition temperature. On the other hand, 
A is a measure of the electron-phonon coupling strength 



A 



dfi 



a 2 F(n) 



(70) 



The behavior of T c with A is very non-linear: For small 
values of A, T c grows exponentially; However, as A in- 
creases, the superconducting transition temperature sat- 
urates. The parameter f2i og is a weighted average of the 
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phonon frequencies 



^io g = exp 



~ j dn log(^)- 



n 



(71) 



and is of the order of the Debye frequency of the material. 
Finally the Eliashberg spectral function is the electron- 
phonon coupling constant averaged on the Fermi surface 



a 2 F(n) = 



N(0) 



EE 

*j A,q 



(72) 

It is widely accepted that McMillan's formula gives a 
quite accurate description of the transition temperature 
for simple, BCS-like, superconductors. Therefore, by 
comparing expressions i)68|) and 169( 1 for the phonon-only 
case, i.e. fi* — 0, we obtain that for BCS-like supercon- 
ductors 



2V(0)£(0) 



-A 



Z(0) « A. 



(73) 



This is an extremely important property of the exchange- 
correlation kernel, which should be fulfilled by any ap- 
proximate functional. 



VII. APPROXIMATIONS TO THE 
ANOMALOUS HARTREE EXCHANGE 
CORRELATION KERNEL 

From the perturbative expansion of the exchange- 
correlation free-energy it is clear that we can split the 
free-energy in three parts: The first contains the purely 
electronic terms, i.e., the terms that do not contain ex- 
plicitly the electron-phonon coupling constant; The sec- 
ond, terms only involving the electron-phonon coupling 
constant; And the last, which we define as the total 
free-energy minus the two first parts, will have mixed 
contributions from the Coulomb and electron-phonon in- 
teractions. The exchange-correlation potentials and the 
exchange-correlation kernels can be split in the same way. 

In this section we develop exchange-correlation kernels 
to be used in the linearized gap equation l|35ll ■ Function- 
al that can be used in the non-linear gap equation H32|) 
are discussed later. This section is divided in two parts: 
First we look at the purely electron-phonon contributions 
to the exchange-correlation kernel. Such functionals are 
developed using the machinery of Kohn-Sham perturba- 
tion theory together with the chain rule introduced ear- 
lier. In the second part, we turn our attention to the 
purely electronic part of the kernel. Two functionals will 
be presented: the first has the form of an local density ap- 
proximation (LDA) , while the second is a functional that 
avoids the direct computation of the screened Coulomb 
matrix elements. The mixed contributions appearing in 
the perturbational expansion of the free energy are ne- 
glected in the current treatment. 
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FIG. 3: The function -jV(0)/C ph (£, f) for Al, Nb and Pb, 
calculated at T = K 



A. Electron-phonon contributions 

In first order in g 2 there are two terms stemming from 
the electron-phonon interaction that contribute to the 
exchange-correlation free energy: given by equa- 

tion {52), andi# given by equation ((56(1 . The exchange- 



correlation kernel derived from F, 
has the form 



(b) 



is non-diagonal and 



K. 



P h 



tanh 



tanh 



( P \ E \9x,q 



To gain further insight into this term, we use a sim- 
plified model: we approximate the electron-phonon cou- 
pling constants by their average value at the Fermi sur- 
face and the electronic energy dispersion is replaced by 
the free-electron model. In Fig. [21 we depict the diagonal 
/C ph (£,£) for aluminum, niobium and lead at zero tem- 
perature for this simplified model. As this contribution 
to the exchange-correlation kernel exhibits particle-hole 
symmetry we only plot the region £ > 0. This term is 
sharply peaked at the Fermi energy (note the logarith- 
mic scale in the £-axis). Furthermore, the width of the 
curves for each material is of the order of the correspond- 
ing Debye frequency. The value of the kernel at the Fermi 
energy can be calculated analytically 



iV(0)/C ph (0, 0) = - [dn a 2 F(Q) 



2 



1 — -^-cotanh 



(75) 



At zero temperature, the value of N(0)IC ph (0, 0) reduces 
to — A, which is the value expected from the comparison 
to McMillan's formula [cf. Eq. 1)73(1] . However, at higher 
temperature 7V(0)/C ph (0, 0) decreases monotonically. 
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The second phononic contribution to the exchange- is a diagonal term, that reads 
correlation kernel coming from the Kohn-Sham pertur- 
bation theory (PT) originates from the diagram i4 c) . it 



rPh.PT _ 



V 13/2 

^3 cosh*(#iy) 



1 



(3/2 



sinh 

1 



(§&) cosh(f^) 



tanh ^ 



9x, c 



2 1 



- n*,,) - /(£.. : ; .«>x (/ : - 2j'te,ei,^A, 9 ) f , (76) 





Al 


Nb 


Mo 


Ta 


V 


Pb 


A^ ph + _Z ph ' sym 


5.59 


15.7 


4.14 


8.48 


23.2 


8.12 


/C ph + Z ph 


7.10 


23.0 


5.23 


11.7 


34.2 


12.8 


Eliashberg 


9.75 


24.7 


7.31 


14.0 


36.4 


12.2 



TABLE I: Transition temperatures from numerical solutions 
of the phonon-only DFT and Eliashberg equations. All tem- 
peratures are in Kelvin. 



where the function I' is defined as 



a 



i'(ti,i j ,n x , q ) = —i(z i ,( j ,n x , q ) 



(77) 



If we try to apply the simplified model presented earlier 
we find that Zf ' PT diverges logarithmically. This di- 
vergence can be traced back to the substitution of <jV 
by its value at the Fermi surface. This problem can be 
solved by retaining the full dependence of the electron- 
phonon coupling constant on the indexes i, and j: 
then decays as a function of energy thereby making the 
integrals present in Eq. I|7t>|l convergent. A closer analy- 
sis of the expressions also reveals that the divergent part 
of the integrands is antisymmetric around the Fermi sur- 
face. Therefore, the divergent integrals would vanish in 
the case of particle-hole symmetry. It seems then reason- 
able to neglect the antisymmetric part of the integrands, 
retaining only the symmetric part. The new functional 
reads 



?ph, sym 



1 



tanh 



x + . (78) 



In Fig. ^ this term is plotted for niobium for several 
temperatures. It turns out that the function Z ph ' sym (£i) 
is a smooth function of the energy, and its value at the 
Fermi surface = 0) is approximately 2A. This is twice 
the value expected from the comparison to McMillan's 
formula [cf. Eq. J73J)]. Furthermore, a careful analy- 
sis of Fig. ^] suggest that Z ph ' sym (£i) can be written as 



is 

a 



'I 1 1 'I 1 1 'I 1 

T = 0.3K 

■,\ T = 0.5K 4 

\\ T = IK - -- 

'■■ x T = 5K 

v • • \ T= 10 K -•- 




le-05 le-04 0.001 0.01 0.1 

Cfe [eV] 



FIG. 4: The dependence of Zf ' sym on temperature for nio- 
bium. 



the sum of two terms: i) one broader and very weekly 
temperature dependent; ii) a second contribution whose 
width decreases significantly with the temperature. Both 
terms contribute with approximately A to the value of 
Z ph ' sym (£i) at the Fermi surface. As the renormaliza- 
tion term Z ph ' sym (£i) appears to be too large, one can 
expect that transition temperatures calculated with this 
functional will be too small. The situation should be 
worst for the strong-coupling superconductors like nio- 
bium or lead, where the renormalization is large. This is 
confirmed by Table^where we list the transition temper- 
atures obtained with the phononic part of the functional. 
These numbers are compared to solutions of Eliashberg's 
equation where we neglected the electron-electron repul- 
sion. 

We believe that the shortcomings of this functional can 
be traced back to the following: Migdal's theorem tells us 
that, to a very good approximation, we can neglect in the 
perturbative expansion diagrams including vertex correc- 
tions due to the electron-phonon interaction. However, 
diagrams including self-energy insertions of phononic ori- 
gin should be included to have a consistent description 
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of the electron-phonon interaction. Therefore the bare 
Green's functions entering in the diagrams depicted in 
Fig. EhiC should be replaced by dressed propagators. In 
a first step to improve our functionals we dressed the 
propagators with a subset of phonon self-energy inser- 
tions. We found that the non-diagonal term /C ph is ba- 
sically insensitive, while the term Z ph ' sym is reduced by 
roughly 20%. This is almost half the correction neces- 
sary to satisfy Eq. JJSJ. We expect that the other 30% 
are accounted for by the remaining self-energy insertions. 
However, this approach is quite involved numerically, so 
we choose a different path to improve our functional. 

We know that the phonon renormalization term should 
have the value A at the Fermi surface. Furthermore, 
this term should have a width comparable to the De- 
bye frequency. It is clear that the broader contribution 
to Z ph ' sym (£i) obeys these requirements. We therefore 
propose to separate the two parts of Z ph ' sym (£j) and use 
the part i) as our renormalization term. We believe that 
this procedure is at least partially justified by the results 
obtained by dressing the Green's functions. The func- 
tional corrected in this way reads 



rPh 



tanh f§£i) j X,q 

[Jfo,Zj,n\, q ) + J(Zi,-Si,tox, q )] , (79) 

where the function J is defined by 

j(t n) = n) - J(£, -n) . (so) 

Finally we have 



J(£,?,Sl) = - 



f (O + n (Q) 



e - e - n 



e - e - n 



(81) 



The functional Z ph is smooth both as a function of the 
energy and as a function of the temperature. Further- 
more, it has approximately the value A at the Fermi sur- 
face. The functional l|79l) . together with the phononic 
kernel J73J, is a central result of our work. It is the func- 
tional that will be used in the calculations of II. In Ta- 
ble[IJ we present the phonon-only transition temperatures 
calculated with this functional. All T c 's are in quite good 
agreement with transition temperatures calculated from 
Eliashberg's equation. We emphasize that the transi- 
tion temperatures in Tableware given exclusively for the 
purpose of testing and/or calibrating the approximations 
made for the phononic part of the exchange-correlation 
functional. T c 's resulting from setting fi* = in the 
Eliashberg equations and setting the Coulomb terms to 
zero in the DFT context have, of course, nothing to do 
with the T c 's observed in nature. For results including 
the Coulomb terms, we refer the reader to II. 



B. Electron-electron contributions 

We now develop functionals that take into account the 
Coulombic part of the interaction. There are two terms 
in the energy functional that give contributions to the 
linearized gap-equation (|35(l . The first is the anomalous 
contribution to the Hartree energy, given Eq. (J2DJ, and 
the second is the exchange term F£ c depicted in Fig. |2Ji. 
The interaction that enters these expressions is the bare 
Coulomb interaction, l/\r — r'\. However, electrons in 
a metal do not feel the bare Coulomb interaction, but 
a much weaker interaction, screened by the sea of elec- 
trons. To take this into account, we take a step back, 
and propose an alternative form to the energy functional 
based on the superconducting version of the local density 
approximation (LDA) 15 . In this approach the exchange- 
correlation energy of the inhomogeneous system is writ- 
ten in terms of the exchange-correlation energy density 
of the homogeneous superconducting electron-gas 



SCLDAr 



i(R), x (R,k)} = 



d 3 R / x h c om Kx(fc)]| „=„ ( r) , (82) 
x=xw(RM) 



where xw(R>k) is the Wigner transform of the anoma- 
lous density of the inhomogeneous system, given by 



Xw{R,k) 



d°se ih8 X [R 



,R 



(83) 



It is easy to see that this definition reduces to the 
usual LDA for non-superconducting systems in the limit 
\ — ► 0. Moreover, it is possible to prove that this is the 
only consistent definition of an LDA for the supercon- 
ducting stated. As an approximation to the exchange- 
correlation energy of the electron-gas one could take the 
random phase approximation (RPA) functional proposed 
in Ref. yj| However, this functional has only been eval- 
uated for a very simple class of pair potentials, namely 
Gaussians centered at the Fermi surface. We therefore 
propose an alternative and simpler form to the Coulom- 
bic contribution to _F XC . For convenience, we approxi- 
mate together the anomalous Hartree and the exchange- 
correlation contributions. Our approximation reads 



/H^W(n)-/£r' NS (n) 



d 3 (r-r')\ X (r-r')\ 2 v TF ( 



(84) 



where v TF (r — r') is the Coulomb interaction screened by 
a Thomas- Fermi model. In coordinate space the Thomas- 
Fermi interaction reads 



v TF (r-r') = 



,-feTF|l — r'\ 



(85) 



with the Thomas-Fermi screening length, &TF, given by 

(86) 



k^ F = 4ttA(0) . 
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By inserting expression (|84|l in the definition of the LDA, 
equation l|82[) . we can identify this approximation as a 
Thomas-Fermi screened anomalous Hartree contribution 
to the free-energy. 

The anomalous Hartree exchange-correlation kernel 
stemming from this term is simply 



K. 



TF 



TF 



(87) 



where the matrix elements of the Thomas-Fermi interac- 
tion are defined by 



TF 



= d 3 r dV <p*(r)<pi(r')v TF (r - r')^(rW*(r') 



(88) 

where the ip^s are the Kohn-Sham orbitals of the inho- 
mogeneous system at hand. 

In II we will compare the results obtained with the 
above approximation with further simplified expressions. 
In the simplest model, the Kohn-Sham orbitals are taken 
to be plane waves with a parabolic dispersion. In this 
case, the kernel can be written in energy space (after 
averaging over the angles) as: 



(k + k'f 



(fc- 



"*TF 



(89) 



with k = y/2(£ - fi) and k' = y/2(£' - /i). Using the 
BCS-like two- well model one can extract the counterpart 
of the Coulomb pseudo-potential /i* from Eliashberg the- 
ory. A crude estimate for r s = 2 gives a value around 0.1, 
which compares well with the typical values of //* for 
simple metals (/i* = 0.10 - 0.15)2£. It should be stressed 
again at this point, that the present method does not 
require fi*. The estimates given here are used to demon- 
strate to which values of [i* our ab initio Coulomb terms 
correspond to. 

While the replacement of the Kohn-Sham orbitals in 
(|88|l by plane waves may be acceptable for simple met- 
als, it will be too crude for more complex materials. In 
those cases it is still possible to avoid the direct com- 
putation of the screened Coulomb matrix elements (|88H . 
by going along the lines described by Sham and Kohn 39 . 
We briefly outline here the main points of this classical 
paper, which deals with an approximate way of getting 
an electron self-energy for the normal state. We assume, 
as usual within the LDA, that our system can be de- 
scribed around the point r by a homogeneous electron 
gas of density n(r). The wave- functions of this electron 
gas can be locally expressed as plane- waves of momentum 
p(r) whose value is determined, in a semi-classical way, 
from the electron energy of the real system. In the sim- 
plest form, the mapping can be obtained from Eqs. (4.5) 
and (4.13) of Ref.|3aas 



El 

2 



& + Hh{n(r)) , 



(90) 



where Hh{n) is the chemical potential of an non- 
interacting homogeneous electron gas with density n. 
Furthermore, we approximate ph(n(r)) by the constant 
Hh(n), where n is the average density of the material. 
We suggest here to approximate the Coulomb interac- 
tion kernel between electrons at energies and £&< by 
the corresponding quantities in the free-electron gas. We 
then replace p 2 /2 — ► & + Hh = Vi> an d rewrite the inter- 
action l|89(l as 



IC SK = 



log 



Vi + r h + 2 VW% 



^tf/ 2 



r\i + rjj - 2 y /rj~rjj + fc TF /2 



In principle, one could consider not onlyp but also fcxF as 
locally dependent on the density n(r). In our simplified 
approach, however, we fix the Thomas-Fermi screening 
length to a constant value. 

Eq. 1|90[) is conceived in terms of wave packets, and is 
valid if n(r) does not vary too much on the scale of the 
Fermi length, exactly as in the normal state LDA. One 
can speculate, however, that when applied to the super- 
conducting state the relevant length scale becomes the 
coherence length, normally much larger than the atomic 
scale. Therefore, we may assume that local variations of 
the density on the atomic scale will not affect the final 
superconducting properties. 

It should be noted that this approximation, although 
derived in the spirit of an LDA, is not a local density 
approximation, since it does not depend explicitly on the 
densities, but implicitly via the single-particle energies £j. 



VIII. FUNCTIONALS FOR THE NON-LINEAR 
GAP EQUATION 



In this section we provide approximations to the 
exchange-correlation kernel that can be used in the non- 
linear gap-equation (|3*2*|l . These functionals will obey to 
one constraint, namely that upon linearization they will 
reduce to the functionals presented in the previous sec- 
tion. This assures that the gap functions obtained from 
Eq. (|32|l and the transition temperatures calculated from 
Eq. (|32|l are consistent. Furthermore, we require these 
functionals to be "well behaved", i.e., without disconti- 
nuities or any other kind of pathological behaviors. 

The simplest way to derive an exchange-correlation 
functional is to use the expressions derived through 
Kohn-Sham perturbation theory in the definition l|18b|) . 

For example, the first phononic contribution 
(Fig. yields the contribution 
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A ph, (b) = 1 z l \ - \p 

^A^J jl A.,, 



2 AjA* + a; A, j e 



y 1 



n 



1 Lf 

yP ?J 



1 x - \ -> Aj 
2^2^ yo^ 

J A,q 





to 




1 



[/(Ej.Bl^A,,)-/^,-^,^,,)] 



f (A,A*-A*A,) | 2 A,^ 2 



tanh cosh 2 (f-Kj) ' Ai Si 

+ (AjA* + A*A,) [I'{E U Ei,tl Kq ) - I'(E t , Q A , g )] 1 . (92) 



It can be seen here that the non-linear gap equation (|32ll , 
in general, does not have the simple structure of a BCS- 
like gap equation and thus goes beyond the simple picture 
of an effective interaction mediating the pairing. How- 
ever, this approach encompasses several problems. First, 
the resulting functionals have extremely complicated an- 
alytical structures and are very hard to interpret in sim- 
ple physical terms. Furthermore, these functionals con- 
tain several divergences and pathological behaviors that 
have to be taken care of. For the time being, we re- 
strict ourselves to using the partially linearized exchange- 
correlation potential, leading to the BCS-type gap equa- 
tion 

i tanh (§Ej) 
A, - - ^ E -^hxc ij Aj . (93) 

3 3 

where J-"hxc i,j are the linearized functionals defined in 
Eq. (|34|) and derived in detail in the previous section. It 
turns out that superconducting gap functions obtained 
with these functionals are in rather good agreement with 
experimental results (see II). 



IX. CONCLUSIONS 

In this work we have developed a truly ab-initio ap- 
proach to superconductivity. No adjustable parame- 
ters appear in the theory. The key feature is that the 
electron-phonon interaction and the Coulombic electron- 
electron repulsion are treated on the same footing. This 
is achieved within a density-functional-type framework. 
Three "densities" : the ordinary electronic density, the su- 
perconducting order parameter, and the diagonal of the 
nuclear TV-body density matrix were identified as suit- 
able quantities to formulate the density-functional frame- 
work. The formalism leads to a set of Kohn-Sham equa- 
tions for the electrons and the nuclei. The electronic 
Kohn-Sham equations have the structure of Bogoliubov- 
de Gennes equations but, in contrast to the latter, they 



incorporate normal and anomalous xc potentials. Like- 
wise, the Kohn-Sham equation describing the nuclear 
motion contains, besides the bare nuclear Coulomb re- 
pulsion, an exchange-correlation interaction. The lat- 
ter is an iV-body interaction, i.e., the nuclear Kohn- 
Sham equation is an iV-body Schrodinger equation. The 
exchange-correlation potentials are functional derivatives 
of a universal functional F xc [n, x, T] which represents the 
exchange-correlation part of the free energy. Approxima- 
tions for this functional were then derived by many-body 
perturbation theory. To this end, the effective nuclear 
interaction was expanded to second order in the displace- 
ments from the nuclear equilibrium positions. By intro- 
ducing the usual collective (phonon) coordinates, the nu- 
clear Kohn-Sham equation is then transformed into a set 
of harmonic oscillator equations describing independent 
phonons. These non-interacting phonons, together with 
non-interacting but superconducting (Kohn-Sham) elec- 
trons serve as unperturbed system for a Gorling-Levy- 
type perturbative expansion of F xc . The electron-phonon 
interaction and the bare electronic Coulomb repulsion, as 
well as some residual exchange-correlation potentials are 
treated as the perturbation. In this way, both Coulom- 
bic and electron-phonon couplings are fully incorporated. 
The solution of the KS-Bogoliubov-de Gennes equation 
(or the KS gap equation together with the normal-state 
Schrodinger equation) fully determines the Kohn-Sham 
system. Therefore, within the usual approximation to 
calculate observables from the Kohn-Sham system, one 
can apply the full variety of expressions for physical 
quantities, known from phenomenological Bogoliubov-de 
Gennes theory, also in the present framework. This ap- 
proach was already successfully applied within a semi- 
phenomenological parametrization of the xc-functional , 
e.g., to the specific hea1*i£ and to the penetration depth^S 
of the cuprates. It should further be emphasized that the 
formalism, developed in this paper, is not restricted to 
perfect periodic systems. It was for this purpose that we 
presented all formulae in terms of general quantum num- 
bers. The formalism can, in principle, be applied as well 
to inhomogeneous systems, containing, e.g., impurities or 
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surfaces, as to perfect periodic crystals. 

In the succeeding paper (II) we will detail the numer- 
ical implementation of this theory and present the first 
full-scale applications to simple metals. 
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Hamiltonian l|29|) for the nuclei gives rise to a set of equi- 
librium coordinates and phonon eigenstates. It can then 
(up to harmonic order) be written as Eq. (jSOJ. The equi- 
librium positions Rq can be employed to define an elec- 
tronic BO Hamiltonian with a lattice potential referring 
to these coordinates 

^BO = f e + (jee + yUtt _ ^3) 

This BO Hamiltonian, without the electron-phonon cou- 
pling, gives rise to the electronic density ur . We now 
add and subtract all Hartree and exchange-correlation 
terms as well as the BO lattice potential to the full Hamil- 
tonian. 

H = (T n + U™ + V£J 

+ (T e + ^ + V^c + A e Hxc ) 

+ L> cc + (L> cn - V£f) 

VScc - 4c - A^c • (A4) 



APPENDIX A: ON THE ELECTRON-PHONON 
COUPLING 

In this appendix we discuss the electron-phonon cou- 
pling potential which appears in the phononic exchange- 
correlation terms. If we decomposed the Hamiltonian 
into the ionic Kohn-Sham Hamiltonian and the electronic 
Kohn-Sham Hamiltonian as the reference system and the 
rest as perturbation, this perturbation would include the 
bare electron-ion interaction. Clearly, the use of the bare 
vertex, i.e. the gradient of the bare nuclear potential 
with respect to the nuclear positions, would yield an 
unphysical electron-phonon interaction. This bare ver- 
tex will be screened by the conduction electrons. This 
screening could be taken care of by a diagrammatical 
resummation2£. 

Here we will sketch a different approach which directly 
generates the screened coupling potential. A natural 
coupling potential in the context of DFT is the gradi- 
ent with respect to the nuclear positions of the effective 
Kohn-Sham potential within the Born-Oppenheimer ap- 
proximation. This is also exactly the quantity which is 
obtained from the standard DFT electron-phonon calcu- 
lations based on linear response theory with respect to 
small lattice distortions. 



Vfl« s ,R(r) 



d 3 ~' 



dV / Hx c(r, r')X(r', r") [v^'(r") 

(Al) 

X(r, r') denotes the linear density-density response func- 
tion and 

/H*c(r,r') = ^-py {vH[n](r)+t^[n](r)} . (A2) 

We are going to outline an approach which generates 
exactly this gradient of the effective Kohn-Sham poten- 
tial as the coupling potential. The effective Kohn-Sham 



The first three terms of this operator represent the nu- 
clear Kohn-Sham Hamiltonian. Assuming that the equi- 
librium density n(r) entering the functional w^ xc [n](r) 
will be close to the equilibrium density (r) result- 
ing from the BO Hamiltonian (|A3|) we can expand the 
Hartree exchange-correlation potential around the BO 
density, 

"HxcMO) = 

>HcJ» + J dV / Hxc W(r,r>n(r') . (A5) 



"Hxcl 



where the small density change Sn(r) is induced by the 
difference of the full electron-ion interaction and the BO 
potential. The density change can, in principle, be calcu- 
lated via linear response to that perturbation. We expect 
this density change to be close to 



6n{r')= d 3 r"X(r',r") V fi <'(r") 



(A6) 



If we keep only the BO part of the electronic Hartree 
exchange-correlation potential, i.e., the term stemming 
from n.R (r), in the electronic Kohn-Sham Hamiltonian, 
we can combine the remainder with the electron-ion in- 
teraction, and can identify (up to first order) 



E 



r Ri 



E 



r R 



0,! 



+ / dV Mx C [nnJ(r,r')5n(r') w VRV SjBO (r). (A7) 

This is the desired result, which allowed us to use the 
electron-phonon couplings, determined by linear response 
calculations, as the coupling potentials in our Kohn- 
Sham perturbation theory. 
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